#include "TROOT.h"
#include "TSystem.h"
#include "TRolke.h"
#include "Riostream.h"
#include "math.h"

void migdalUncertainty()
{
    //
// Gaussian uncertainty in the background estimate
// Gaussian  uncertainty in the efficiency estimate
//
    cout << endl<<"======================================================== " <<endl;
    Double_t bm = 0.25;     // expected background
    Int_t x = 6;          // events in the signal region
    Double_t sdb = 0.0724;    // standard deviation in background estimate (stat+sys)
    Double_t em = 0.144;     //  measured efficiency
    Double_t sde = 0.0190263;    // standard deviation of efficiency (stat+sys)

    // make TRolke objects
    TRolke tr;   //

    Double_t ul ; // upper limit
    Double_t ll ; // lower limit

    tr.SetCLSigmas(1);

    tr.SetGaussBkgGaussEff(x,bm,em,sde,sdb);
    tr.GetLimits(ll,ul);
    cout << "For model 3 : Gaussian / Gaussian" << endl;
    cout << "the Profile Likelihood interval is :" << endl;
    cout << "[" << ll << "," << ul << "]" << endl;

    cout << "***************************************" << endl;

    Int_t num_NR = 8.17e5;

    Double_t center_val = (x - bm) / Double_t(em * num_NR);
    Double_t lo_err = -ll / num_NR + center_val;
    Double_t hi_err =  ul / num_NR - center_val;
    
    Double_t syserr_NR = 40101.;
    Double_t err_NR = sqrt(num_NR + syserr_NR*syserr_NR);

    lo_err = sqrt( pow( lo_err/center_val, 2) + pow(err_NR/num_NR, 2)) * center_val;
    hi_err = sqrt( pow( hi_err/center_val, 2) + pow(err_NR/num_NR, 2)) * center_val;

    cout << "Center value and Asymmetrical uncertainties (Including both stat and sys err):" << endl;
    cout << "( " << center_val*1e5 << " - " << lo_err*1e5 << " + " << hi_err*1e5 << " ) * 10^{-5}" << endl;
    cout << "***************************************" << endl;

    cout << endl<<"======================================================== " <<endl;

    cout << "To calcualte the siginifance, we don't need to consider the stderr of efficiency." << endl;
    cout << "Thus we use the model 5." << endl;
    cout << "For model 5 : Gaussian / Known" << endl;
    Double_t e = em;    // efficiency
    tr.SetCLSigmas(1);
    tr.SetGaussBkgKnownEff(x,bm,sdb,e);
    tr.GetLimits(ll,ul);
    //cout <<  "the Profile Likelihood interval is :" << endl;
    //cout << "[" << ll << "," << ul << "]" << endl;
    cout << "***************************************" << endl;
    for(Int_t sig=1; sig < 8; ++sig)
    {
        tr.SetCLSigmas(sig);
        Int_t ncrt;
        tr.GetCriticalNumber(ncrt);
        cout << "Critical number for " << sig << " sigma: " << ncrt << endl;
    }

}
